Gravitational Waves from Nonlinear Couplings of Radial and Polar Nonradial Modes 

in Relativistic Stars 
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The post-bounce oscillations of newly-born relativistic stars are expected to lead to gravitational- 
wave emission through the excitation of nonradial oscillation modes. At the same time, the star 
is oscillating in its radial modes, with a central density variation that can reach several percent. 
Nonlinear couplings between radial oscillations and polar nonradial modes lead to the appearance 
of combination frequencies (sums and differences of the linear mode frequencies). We study such 
combination frequencies using a gauge- invariant perturbative formalism, which includes bilinear 
coupling terms between different oscillation modes. For typical values of the energy stored in each 
mode we find that gravitational waves emitted at combination frequencies could become detectable 
in galactic core-collapse supernovae with advanced interferometric or wide-band resonant detectors. 

PACS numbers: 04.30.Db, 04.40.Dg, 95.30.Sf, 97.10.Sj 



I. INTRODUCTION 

Oscillations of neutron stars can arise in strongly non- 
linear phases of stellar evolution, such as during the post- 
bounce phase in core collapse of massive stars or during 
the merger of compact binary systems. Linear pertur- 
bation theory is appropriate for describing the dominant 
properties of stellar pulsations, but neglects several in- 
teresting nonlinear effects that can modify or enrich the 
linear description. For instance, a nonlinear approach is 
required for understanding how energy transfer between 
various classes of modes could saturate the r- and f-mode 
instabilities in rotating stellar models (see [l| and refer- 
ences therein), or limit the persistence of the bar- mode 
instability, through nonlinear coupling between m — 1 
and m = 2 modes 0, H, Nonlinear mode couplings 
also become important in the development of a g-mode 
instability, that was recently observed in core-collapse 
simulations [s'] . Another nonlinear effect relevant for pul- 
sations in rapidly rotating stars is the mass-shedding- 
induced damping of pulsations [H, 01 • 

The main nonlinear effect appearing even for mildly 
nonlinear pulsations is the presence of combination tones, 
which are nonlinear harmonics whose frequency appears 
as a linear combination of linear eigenfrequencies (see 
0) H) ) ■ The number of linear terms in the combination 
tones is equal to the nonlinear order of the mode cou- 
pling. In perturbation theory, these nonlinear features 
can be addressed starting at bilinear order, for the cou- 
pling of two linear modes. 

The presence of nonlinear harmonics can lead to ex- 
change of energy between modes or even to two-mode or 
three-mode resonances or parametric instabilities Q . In 
rotating proto- neutron stars, rotational effects can in- 
crease the number of possible resonances between ax- 
isymmetric modes Q. The parameter space for pos- 
sible resonances or instabilities is even larger for non- 



axisymmetric oscillations. In this case, each Z-mode is 
split in its m components proportionally to the rotation 
rate of the star. Therefore, at large rotation rates, several 
modes could satisfy resonance conditions • 

Gravitational waves emitted by neutron star pulsations 
are interesting sources for current detectors, but the high- 
frequency band will only start becoming accessible with 
the construction of advanced detectors, such as Advanced 
LIGO 11], VIRGO+ and Advanced VIRGO [H, HI. 
There are also plans for detectors with significantly in- 
creased sensitivity at several kHz, such as the GEO-HF 
project [Til and the proposal for the Dual wide-band bar 
detector [15[ . Nonlinear effects in neutron star pulsations 
could be probed with such instruments. 

Interest in the nonlinear dynamics of stellar pulsations 
exists also for variable stars, such as Cepheids, RR Lyrae 
and S Scuti stars. The modulation of the pulsation am- 
plitude (Blazkho effect) or other irregular oscillations, 
which have been observed in the velocity and light curve 
of these stars, could be explained by nonlinear interaction 
between various modes. These studies have been carried 
out mainly in Newtonian perturbation theory by using 
amplitude equations (see e.g. [H, [l3, Ull and references 
therein) . 

General Relativity provides a more appropriate frame- 
work for studying compact objects, where relativistic ef- 
fects, such as the dragging of the inertial frames or the dy- 
namical role of the spacetime, can influence the spectral 
properties of hydrodynamical modes or introduce new 
features, such as w- modes [l^. Modern high-resolution 
methods [1^ (in particular, the 3rd-order PPM scheme) 
implemented in numerical relativity codes, have allowed 
in recent years the study of nonlinear pulsations in non- 
rotating and rotating stars (see 0, 0, (H [H HI ) . The 
above studies included 2-D and 3-D simulations of ax- 
isymmetric modes in either the relativistic Cowling ap- 
proximation, the spatially conformal flatness (CFG) ap- 
proximation or full general relativity. Nonlinear rela- 
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tivistic radial pulsations have also been addressed with 
a different method in [l^, , where the nonlinear dy- 
namics was studied as a deviation from an equilibrium 
background. 

Complementary to fully nonlinear approaches, nonlin- 
ear perturbative techniques can be very useful in weakly 
and mildly nonlinear regimes for understanding how a 
physical process is affected by nonlinear effects. In ad- 
dition, the dimensional reduction of the problem by sep- 
aration of variables (when possible) yields a much less 
computationally expensive approach, suitable for large 
parameter studies. In order to correctly address the 
gauge freedom of general relativity, a relativistic pertur- 
bation theory has been introduced in [26| . and then ex- 
tended to physical systems where two or more variables 
can be assumed as perturbative parameters [13, HE HI] ■ 
For spherically symmetric and time dependent space- 
times, a gauge invariant formalism for studying spher- 
ical and nonspherical perturbations was introduced by 
Gerlach and Sengupta l30l . [sij . This formalism was fur- 
ther developed in [32, l33| and then extended (with par- 
tial gauge-invariance at second and higher perturbative 
orders) in [sj] (from now on GSGM formalism). Efforts 
are now under way for formulating a fully gauge- invariant 
theory at second or higher order [34]. 

Second-order perturbative studies have been car- 
ried out in cosmology [IH] , or for investigating the 
Schwarzschild 36., .37, HI] and Kerr [s^ spacetimes. In 
the two-parameter relativistic perturbative framework, 
the coupling between the radial and nonradial oscillations 
of perfect-fluid, spherically symmetric neutron stars has 
been recently investigated in ll,|40,|4l[. In particular, the 
coupling between the radial and axial oscillations exhibits 
an interesting resonance when a radial mode frequency 
approaches the frequency of an axial w-mode 

In this paper, we address the coupling between the ra- 
dial and polar nonradial modes by numerically solving 
the perturbation equations introduced in [4^ for poly- 
tropic, nonrotating relativistic stars. In the polar sector, 
new interesting results are expected as the spectrum of 
polar modes is generally richer than the axial case. In 
particular, we set up initial conditions (in following recent 
numerical simulations of core collapse [il,|43[) for study- 
ing pulsating protoneutron stars. In the post-bounce 
phase a newly-born neutron star is expected to mainly 
oscillate in its radial (I = 0) and nonradial quadrupole 
{I = 2) modes [IE 113 with a variation in its central 
density of several percent. After investigating several 
initial pulsating configurations, we have found that for 
a 5% variation in the central energy density and about 
lO^^M© energy stored in nonradial oscillations some bi- 
linear combination tones are within the sensitivity win- 
dow of advanced and newly proposed gravitational wave 
detectors. A possible detection of these combination 
tones, in addition to the linear modes, would provide new 
important information for solving the high-density equa- 
tion of state puzzle through gravitational-wave astero- 
seismology. Note that the perturbative equations as de- 



veloped in [1, [13, nil do not contain back-reaction terms. 
This issues is the subject of future work. 

The plan of this paper is as follows. Sections [TTl and Hill 
are dedicated to the description of the background stel- 
lar configuration and the perturbative framework, respec- 
tively. The set of perturbed initial data, boundary con- 
ditions and main properties of the numerical code are 
addressed in Sec. |lVl In Sec. |Vl we present the temporal 
and spectral properties of the fluid and Zerilli variables, 
which arise from the coupling between radial and polar 
nonradial oscillations. The detectability of the gravita- 
tional wave signal is discussed in Sec. lVIl In Appendix El 
we present the equations that complete the system of 
nonradial perturbation equations given in poj . while in 
Appendix [B] we derive the pulsational kinetic energy of 
nonradial oscillations in terms of GSGM formalism. 

In this paper the geometrical units are adopted, where 
G = c= 1. 



II. BACKGROUND EQUILIBRIUM 

The background spacetime is the equilibrium configu- 
ration of a perfect-fluid spherically symmetric nonrotat- 
ing star, which is described by the following spacetime 
metric: 

ds^ = -e^'^dt^ + e^^'dr^ + r^{de^ + sin^ Odc^^) , (1) 

where $ = <I>(r) and A = A(r) are obtained by solv- 
ing the Tolman-Oppenheimer-Volkov (TOV) equations. 
This system of equations is closed by an equation of state 
(EoS) that characterizes the fluid properties. In this pa- 
per we consider a relativistic barotropic EoS: 

p = Kp^, (2) 
£ = P+^, (3) 

where p is the pressure, p and e are the rest mass den- 
sity and the total energy density respectively, and K 
and r are the polytropic parameters. For a central den- 
sity pc = 5.87 X 10~^ km~^ and polytropic parameters 
K = 217.858 km^ and F = 2, the TOV equations pro- 
vide a stellar model with typical mass M = I.4OM0 and 
radius R = 14.151 km 0. 

III. PERTURBATIVE FRAMEWORK 

The system of perturbation equations for studying the 
coupling between the radial and nonradial oscillations of 
spherical star has been introduced in 110, EH . There- 
fore in this paper we only report the main properties of 
the perturbation framework and equations. For more de- 
tails see the references cited above. 

The study of linear perturbations of a spherically 
symmetric background can be significantly simplified by 
adopting the expansion of the perturbed quantities in 
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scalar, vector and tensor harmonics. This technique en- 
ables us to separate in the perturbed quantities the an- 
gular dependence from the time and radial parts. This 
reduces the problem to just one spatial dimension. Any 
harmonic component of the linear perturbative variables 
is then identified by the harmonic indices {£, m) and is dy- 
namically independent. A further sub-classification can 
be carried out on a spherically symmetric background, 
where two perturbation classes with opposite parity and 
independent dynamics can be defined, namely the axial 
[odd-parity) and the polar {even-parity) perturbations. 

The above perturbative technique can be extended to 
second order perturbations. However, due the hierar- 
chical structure of the perturbation methods the non- 
linear perturbations obey inhomogeneous perturbation 
equations, where the homogeneous part has the same dif- 
ferential structure as the linear perturbation equations, 
whereas the source terms are made by the product of 
linear perturbations. From the form of these quadratic 
terms, we can notice that the [i, m) second order per- 
turbations depends on self coupling terms, which couple 
linear perturbations with the same harmonic index and 
mixed terms, which are instead the product of linear per- 
turbations with different harmonic indices. 

In this paper as well as in [§, [33, , we select the cou- 
pling between the radial (£ = 0) and nonradial {I > 2) 
perturbations, and concentrate on the numerical sim- 
ulations of the polar perturbative class. For the cou- 
pling between the radial and axial nonradial oscillations 
see [i, Hlj . In order to study the gauge properties of the 
linear and nonlinear perturbations, we have found con- 
venient to use the 2-parameter relativistic perturbation 
theory [2^ and the gauge invariant formalism of Gerlach 
and Sengupta [sol . Isil [H, , in the version developed 
by [13, , that we call GSGM formalism. We have la- 
belled the radial and nonradial perturbation with two 
distinct parameters and studied the gauge invariance of 
our nonlinear quantities. This particular coupling can 
then be studied by means of variables that are gauge 
invariant when the linear radial perturbation gauge is 
fixed. Adopting the so-called radial gauge for the radial 
perturbations, the problem becomes essentially that of 
linear polar perturbations on a time-dependent spherical 
background. 

The remainder of this section is dedicated to the de- 
scription of the main properties of the radial and non- 
radial perturbations and their coupling terms. We label 
linear radial quantities by a ^^'^-^ superscript, linear non- 
radial quantities by a ^^'^^ superscript and coupling terms 
by a (1^1) superscript. The radial pulsations are described 
by a set of four perturbative fields, namely two metric 
quantities S'^^''^\ rj^^'^') and two fluid variables, which 
are the enthalpy H^^'^^ and the velocity 7(1'°) variables. 
They obey three first order in time evolution equations 
and two constraints, as there is a single radial degree of 
freedom. We can then solve a hyperbolic-elliptic system 
of equations, which is formed by two hyperbolic equa- 
tions, for --/(I'O) and H'-^'^\ and the Hamiltonian con- 



straint that at any time step updates the variable 5^^'"^ 
while ry^i'*^) is determined by the second elliptic equation. 

The linear polar nonradial oscillations can be studied 
with a set of gauge invariant quantities [H, HI, HO] , the 
three metric perturbations S^°'^\ k^^'^\ i/j^"'^) and the 
three fluid variables ij'"'-'^), -y^'^'i), a^^^'^^ Among the dif- 
ferent systems of perturbation equations that are avail- 
able in the literature, we choose a system of three par- 
tial differential equations for the three variables S^'^'^\ 
/j(o,i) a,nd H'-^-'^\ It consists of two hyperbolic equations, 
which describe respectively the gravitational waves and 
sound waves propagations, and an elliptic equation, i.e 
the Hamiltonian constraint 0, 13, . Since the Hamil- 
tonian constraint is used for updating at any time step 
one of the unknowns of the problem, the errors associ- 
ated with the violation of this constraint are corrected. 
The other nonradial variables ■0^'^'^^ ry(o,i) ^ which 
are necessary for updating the source terms of the second 
order perturbative equations, can be obtained by three 
hyperbolic partial differential equations that are reported 
in Appendix [XI 

The nonlinear coupling terms obey the same system of 
equations as the linear nonradial terms, but with non- 
vanishing source terms in the region interior to the star. 
Therefore, we evolve the two metric quantities S^^'^'^ and 
/c*^^'^^ and the fluid variable H^^'^\ In order to under- 
stand the structure of the perturbative framework we 
define the following set of linear and nonlinear pertur- 
bations: 



^(1,0) ^ (5(1,0)^ ^(1,0)^ ^(1,0)^ ^(i,o)j ^ 

^(0.1) = (5(o.i),fc(04)^^(o,i)^ij(o,i)^^(o,i)^^(o,i)) ^ 



Qii.i) ^ (^5(1,1)^^(1,1)^^(1,1)) ^ 



(5) 
(6) 



and write the system of the three perturbative equations 
for the coupling as follows: 



g(i,i) =s 7^(l'°)®AA(0'l) 



(7) 



where L^n is an operator representing the homogeneous 
part of the system, which has the same differential struc- 
ture as the linear nonradial perturbation equations for 
the variables S^^''^\ k^^'^^ and i/*^^'^^. The source terms 
are instead represented by the operator S, which also 
contains spatial and first order time derivatives. 

Outside the star, the linear polar nonradial perturba- 
tions have just one degree of freedom, which is described 
by means of the Zerilli-Moncrief function [46l|. In terms 
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of the gauge-invariant variables it is is given by ^ 



^(0.1) 



£{£+!) [{£{£ + 1) - 2)r + 6M] 



, (8) 



where M is the mass of the star. The radial-nonradial 
couphng terms have the same angular dependence and 
satisfy an equivalent Zerilli-Moncrief wave-like equation 
as the linear nonradial perturbation. We denote this vari- 
able as Z^^'^); it is constructed from an expression of the 
form of Eq. ([8]) with all the variables of type replaced 
by the corresponding '^'^^ ones. For any multipole {£, m), 
the total Zerilli-Moncrief function is obtained as 



(9) 



where we have explicitly reintroduced the multipolar in- 
dices. As a result, the total emitted power in gravita- 
tional radiation (at infinity) is computed as (see for ex- 
ample \4l\ ) 
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(10) 



where the overdot denotes time derivative. 



IV. SETTING UP NUMERICAL SIMULATIONS 

The numerical code for studying the dynamical evolu- 
tion of polar coupling perturbations has the same struc- 
ture as the code developed in [9, 41] for the axial cou- 
pling case. We had to replace the parts of the code 
that compute the axial linear and coupling perturba- 
tions with the new routines that are appropriate for the 
polar sector. For details about the overall structure of 
the code see |4l| , while for the part that solves the 
polar nonradial perturbations see [3, |4^. This latter 
part has also been used for treating the equations for the 
radial-nonradial terms, by adding the source terms found 
in HIH. 



A. Initial Data 

The independence of linear perturbations from the har- 
monic index £ implies that we need to separately excite 
the radial and nonradial oscillations. In order to sim- 
plify the study of numerical simulations and identify the 



^ In the reference there is a typo in the definition of the Zerilli- 
Moncrief function that here we have corrected. In the same pa- 
per, equation (79) must be corrected as follows: U 
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TABLE L Eigenfrequencies (in kHz) of the first eight fluid 
modes of the radial and nonradial £ = 2 polar perturbations. 
The second column corresponds to the results from calculation 
in the frequency domain (the solution of the Sturm-Liouville 
problem [y]), while the third column displays the results of 
an FFT of the enthalpy time evolution. In square brackets 
we denote the relative difference between the frequencies ex- 
tracted from the time evolution and the Sturm-Liouville solu- 
tions. The fourth column shows the frequencies of the polar 
nonradial oscillations, which are determined by an FFT of a 
time-domain simulation. 



Radial 


Freq. code 


FFT 


Nonradial 


FFT 


F 


1.443 


1.437 [0.4%] 


H 


1.587 


Hi 


3.955 


3.951 [0.1%] 




3.757 


H2 


5.916 


5.913 [0.05%] 


'P2 


5.699 


Hs 


7.775 


7.771 [0.05%] 


^P3 


7.614 


H4 


9.590 


9.587 [0.03%] 


'P4 


9.419 


Hs 


11.38 


11.368 [0.1%] 


'P5 


11.25 


He 


13.15 


13.152 [0.01%] 


^P6 


13.03 


Ht 


14.92 


14.916 [0.02%] 


'P7 


14.781 



y(i,0) 



correct nonlinear harmonics we have excited the linear 
perturbations by selecting specific modes. For radial pul- 
sations we have set up a Sturm-Liouville problem for the 
radial velocity perturbations and solved it numerically 
with relaxation methods. As tested in [1, |4l[, this code 
determines the eigenfrequencies of radial modes with an 
accuracy to better than 0.2% with respect to the pub- 
lished values. With an appropriate choice of the initial 
phases for the various perturbed variables, we can ex- 
cite the radial eigenmodes by providing only the eigen- 
functions associated with the radial velocity perturba- 
tion --/(I'O). The simulations are numerically stable over a 
large multiple of the largest oscillation period that we are 
interested in. Using Fast Fourier Transformations (FFT) 
of the time evolution of selected variables, we can repro- 
duce the linear mode frequencies of radial modes with an 
accuracy to better than 0.4%, as shown in Table IH More 
details can be found in [§, |41| . 

There are various sets of initial conditions for excit- 
ing the polar nonradial perturbations which satisfy the 
Hamiltonian and momentum constraints on the initial 
slice dl, m. Hi]. After having explored all the differ- 
ent choices mentioned in the literature, we have not no- 
ticed any important qualitative difference with respect 
to the coupling perturbations. Therefore, in this paper 
we consider only a representative case, which consists in 
perturbing the fluid with an enthalpy eigenfunctions, im- 
posing S'f^d) = (conformally flat initial data) and de- 
termining A:*-*^'^' by solving the Hamiltonian constraint. 

The enthalpy eigenfunctions are obtained with the 
eigenfunction recycling method developed in @, 0] , where 
the eigenfunction is extracted after a simulation by means 
of FFT transformations at every grid point. In a first 
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FIG. 1: Normalized eigenfunctions of the perturbed enthalpy 
for the fundamental nonradial mode ^f and its first four over- 
tones, ^pi to ^P4 {upper panel). The lower panel displays 
the Fourier power spectral density of the evolution of the 
perturbed enthalpy, in simulations where each of the above 
eigenfunctions was used as an initial perturbation. In all 
cases, the desired mode is predominantly excited, while other 
modes (still excited due to truncation errors) have an ampli- 
tude which is several orders of magnitude smaller. 



simulation, we use trial eigenfunctions that excite var- 
ious modes. The eigenfunction of the desired mode is 
then extracted and used as an initial perturbation for 
a second simulation, in which case the desired mode is 
predominantly excited. This recycling procedure can be 
repeated until the amplitude of the other modes becomes 
sufficiently small. Examples of extracted eigenfunctions 
are shown in the top panel of Fig. [1] which displays the 
eigenfunctions of the fundamental quadrupole mode ^f 
and its first four overtones ^pi to ^p4, for the equilib- 
rium model mentioned in Sec. |TT1 The lower panel dis- 
plays the Fourier power spectral density (PSD) of the 
time evolution of the perturbed enthalpy, in simulations 
where each of the above eigenfunctions was used as an 
initial perturbation. In all cases, we observe that the de- 
sired mode is predominantly excited, while other modes 
(still excited due to truncation errors) have an amplitude 
which is several orders of magnitude smaller. 



B. Case Studies 

The nonlinear coupling between radial and nonradial 
modes is studied here in a number of particular cases, 
in which the perturbed initial data consists of a linear 



TABLE II: Initial configurations (CASE I - IV), constructed 
by various linear combinations of eigenfunctions correspond- 
ing to particular radial and nonradial modes. For the relative 
amplitudes used in these linear combination, see the main 
text. 



CASE 


Radial modes 


Nonradial modes 


I 


F 




II 


F, Hi, H2, H3, H4, H5 


2f 


III 


F 


2f 2„ 2^ 2„ 2^ 2^ 
I, Pi, P2, P3, P4, P5 


IV 


F, Hi, H2, H3, H4, H5 


2f 2„ 2„ 2„ 2„ 2„ 
I, Pi, P2, P3, P4, P5 



combination of various radial and (quadrupole) nonradial 
modes. Table |TT] provides a summary of the modes used 
in each case. In CASE I, only the fundamental radial 
and the fundamental nonradial modes are present in the 
perturbed initial data. CASE II includes, in addition, 
the first five radial overtones. CASE III includes the 
fundamental radial and nonradial modes plus the first 
five nonradial overtones. Finally, CASE IV includes all 
radial and nonradial modes up to the fifth overtone. 

The amplitude of each linear mode is determined by 
choosing a particular value of its pulsational kinetic en- 
ergy. For radial pulsations, the kinetic energy is com- 
puted at an instant of vanishing Lagrangian displace- 
ment (vanishing potential energy, maximum kinetic en- 
ergy) as [sot : 

41^") = 2^ dr (£ + p) 6^^+* (r 7(1'°)) ' , (11) 

where 7(1'°) is related to the perturbed radial velocity. 

For nonradial pulsations, the kinetic energy can be de- 
termined in terms of the GSGM quantities in the Regge- 
Wheeler gauge as follows (see Appendix |B]) : 

i^r^ ^lj\ris+p)e''^^ 

'r^{,i^^^^-^)\iii+l){aiO:^^y , 

(12) 

where 7('^'^) and a'"'^'' are proportional to the radial 
and longitudinal components of the velocity, respectively, 
while the metric perturbation i/j^^'^) is proportional to 
6gro- In order to determine the amplitudes of the ini- 
tial data, we average the pulsational kinetic energy of 
Eqs. pT|) - (fT^ over several oscillation periods. 

After having determined the required amplitude for 
each mode to get a certain pulsation energy, the radial 
pulsations are excited by a linear combination of the vari- 
ables 7!^'°^ (related to the radial velocity eigenfunction), 
where the index n denotes different modes: 

^(1,0)^^^(1,0). (13) 

n 
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FIG. 2: Time evolution {upper panel) and the spectral prop- 
erties {lower panel) of the enthalpy for CASE I, when the 
star is oscillating at first perturbative order in the radial and 
nonradial fundamental modes. In the {upper panel), the en- 
thalpy is given as log |//'''-'^ |, where i,j = 0,1. The radial 
and nonradial enthalpy are shown with solid and 
dashed lines respectively, while H^^'^^ is represented with a 
dotted line. The nonlinear harmonics have been denoted as 
f± = 2f ± F and p* = ^Pn ± F. 



Correspondingly, for nonradial oscillations 
combination of the enthalpy eigenfunctions H 
to initiate the pulsations, is given by 



the 

(0,1) 



i7(o,i) 



linear 
, used 



(14) 



For simplicity we consider vanishing initial data for all 
coupling variables, i.e. Q'^'^-' = 0. The Hamiltonian and 
momentum constraints are then violated on the initial 
Cauchy hypersurface and the numerical solutions will be 
effected by an initial transient of short duration [§, [4l[ . 
This initial transient is not affecting the main results pre- 
sented here. 
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FIG. 3: For the same initial pulsating configuration of Fig. [2] 
(CASE I) the upper panel depicts the total enthalpy H as 
given by Eq. (|15|) in the equatorial plane {6 — tt/2, solid 
line) and in the polar direction {9 = 0, dashed line). Lower 
panel: the Fourier power spectral density of H in the equato- 
rial plane. 



been also implement for polar modes. During the evo- 
lution, some of the junction conditions for the coupling 
variables could not be correctly imposed, due to the mo- 
tion of the perturbed stellar surface (see 0, [131 )• We 
address this issue as follows: i) at the stellar surface, we 
determine the maximal Lagrangian displacement of the 
linear oscillations and ii) we impose the junction condi- 
tions for the coupling variables on a hypersurface which is 
always inside the oscillating star. This method is equiv- 
alent to removing a thin outer layer close to the stel- 
lar surface, with a thickness depending on the pulsation 
amplitude. In practice, for small-amplitude pulsations, 
only very few grid points near the stellar surface of the 
equilibrium model need to be neglected for the coupling 
variables. This does not have a noticeable quantitative 
effect on our main results. 



C. Boundary Conditions 

To complete the description of the initial-value prob- 
lem, one needs to impose the boundary conditions at the 
origin, the stellar surface and at infinity. This topic has 
been addressed in detail in |4l| and hence we do 

not repeat it here. However, it is worth mentioning that 
the surface treatment used in the axial case [9|, |41|] has 



V. RADIAL-NONRADIAL MODE COUPLINGS 

We discuss now the numerical results for the coupling 
of radial and nonradial pulsations, for the specific cases of 
initial configurations described in Sec. IIVBI In order to 
interpret the numerical simulations, we first describe the 
main properties that we expect from the structure of the 
equations for the coupling perturbative terms Q*-^'^' . The 
linear radial oscillations will contain the eigenfrequen- 
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FIG. 4: Gauge-invariant waveforms {top panel) and Fourier 
power spectral density {bottom panel) for CASE I: we show 
together the coupUng Zerilh-Moncrief function Z^^'^^ {dashed 
line) and the total one Z = Z^"'^' + {solid line). 



FIG. 5: Fourier power spectral density of Z^^'^^ {dashed line) 
and of Z {solid line). The top panel refers to CASE II and 
the lower panel to CASE III (see text for details). In order 
to simplify the labelling of the figure we have defined = 
Hn ± ^f , f± = ^f ± F and p± = ^pn ± F. 



cies t'j ' of the particular radial modes excited by the 
initial data. Similarly, the linear nonradial oscillations, 
will contain the frequencies j^I'^'^'' of the particular nonra- 
dial modes. The variables describing the radial-nonradial 
coupling terms, Q*^^'^-*, which obey equations will 
then be driven by two different kinds of frequencies: i) 

the natural frequencies h'n'^\ associated with the homo- 
geneous part of the nonradial equation ([7]), which are now 
excited by the forcing term S, and ii) the nonlinear cou- 
pling frequencies, the so-called combination tones, whose 
frequencies are linear sums or differences of linear mode 
frequencies. For the second-order coupling between ra- 
dial and nonradial oscillations, an example of such a com- 
bination frequency is z^t^^^)^ = jy(i.o) j/(o,i)^ 



A. Case I 

In CASE I, only the fundamental radial and the fun- 
damental £ — 2 nonradial modes are present in the per- 
turbed initial configuration. In order to set the amplitude 
of each mode, we take into account recent numerical sim- 
ulations (sec [42, 43, 51, 52, 53] and references therein) 
which suggest that over a total simulation time of about 
20 ms the total energy emitted in gravitational waves 
can be of order 10~* M©. Furthermore, the post-bounce 
radial pulsations in the above simulations induced a vari- 
ation in the central energy density of the order of 1-5%. 



In our numerical simulations, we find that a gravita- 
tional wave energy of 1.11 x 10"^ Mq (as determined by 
time integration of Eq. (|10p ) is emitted after 20 ms of 
continuous monochromatic emission of the fundamental 
£ = 2 mode, when the latter has an average kinetic en- 
ergy of (Sfc^'^') = 5.0 X 10^^ Mq. In practice, we stop 
our simulation at 34 ms, so that the energy emitted in 
gravitational waves is 1.74 x 10~^ Mq. The radial pul- 
sations are excited with an amplitude that corresponds 
to an average kinetic energy {E^}'^^) — 5.0 x 10"'' Mq, 
which leads to a variation in the central energy density 
of 1%. 

The upper panel of Fig.[2]displays the time evolution of 
the enthalpy variables H''^'°\ and i/^^-i' averaged 

at any time step along the spatial coordinate. While the 
two linear variables are nearly perfectly monochromatic, 
the bilinear term /f^^'^^ contains several frequencies. In 
the PSD of H^-^'^\ shown in the lower panel of Fig. ^ 
several combination frequencies, such as f* = ^f ± F 
and Pjf = ^Pn ± F, for n = 1, 2, 3, can be clearly iden- 
tified. We have also verified that the amplitude of i/*^^'^^ 
is nearly the product of the two oscillation amplitudes 
of 77(1^°) and H'^^'^\ The top panel of Fig. [3] shows the 
corresponding time evolution of the (averaged) total en- 
thalpy H {a± 9 — 7r/2, solid line, and = 0, dashed line) 
defined as 

H ^ + ~ ^ + g^^-^)) , (15) 



8 



where the normaUzation constant of the spherical har- 
monic is contained in the radial variable H'-^-'^\ 
One can notice a modulation of the amplitude, which is 
mainly due to the different angular dependence of the ra- 
dial and nonradial oscillations. The lower panel of Fig. [3] 
shows the PSD oi H at 9 = tt/2. Both linear radial 
and nonradial fundamental modes are present. Of the 
several combination frequencies present in H^^'^\ only 
f ^ = ^f ± F have a large enough amplitude to be clearly 
visible here. We recall that the two lowest-order lin- 
ear pressure mode frequencies seen in the lower panels 
of Fig. [5] and Fig. [3] are present because the initial data 
for the linear fundamental radial mode are not perfect, 
but also contain additional frequencies with very small 
amplitudes. The top panel of Fig. [4] shows the time evo- 
lution of Z^^'^^ (dashed line) and Z (solid line) extracted 
at r = 100 km. We note that the bilinear function Z^^'^'^ 
is about two order of magnitude smaller than Z. In the 
lower panel of the same figure, the corresponding PSD 
are displayed. The main combination frequencies are 
f^ = ^f ± F, whose detectability will be addressed in 
Sec. EI below. 



B. Cases II and III 

Next, we study two initially pulsating configurations 
where several modes are contemporarily excited. In 
CASE II, we excite the nonradial fundamental mode ^f 
and a linear combination of radial modes up to the fifth 
overtone. In CASE III, we excite the radial fundamental 
mode F and a up to the fifth pressure nonradial mode. 
The pulsation amplitudes of the two fundamental modes 
F and ^f have the same value as in CASE I. For the over- 
tones, we assume that the kinetic energy decreases pro- 
portionally with the order of the mode. This assumption 
is motivated by the spectral properties of the gravita- 
tional waves emitted in core collapse simulation. In par- 
ticular, we demand that each overtone or pressure mode 
of order n has one fifth of the energy stored in the n — 1 
mode; i.e., 



(1,0)\ 



{Er 



iiEr 

(E, 



k /Hn- 

(0.1)\^ 



(16) 
(17) 



where for n = the F and ^f modes are implied. Notice 
that due to phase-cancellations, the total energy, when 
exciting several modes, is not a linear sum of the in- 
dividual energies. Using Eq. ([TT|) and Eq. p2|) for the 
total perturbation, we determined that when the radial 
or nonradial oscillations are excited up to the fifth over- 
tone, the radial and nonradial average kinetic energies are 
{Ei^-°^) = 6.25 X 10-^ Mq and (£;^°'^^) = 6.0 x lO'^ M© 
respectively. 

For CASE II, the upper panel of Fig. O displays the 
PSD of the bilinear Zerilli-Moncrief fimction Z'^'^^ as 



well as that of the total Z. On the other hand, the bot- 
tom panel of the same figures exhibits the same informa- 
tion for CASE III. The spectra clearly show the presence 
of several expected combination frequencies: in particu- 
lar, for CASE III, a repetitive triplet structure (such as, 
e.g. p|',p2,P;^) is evident in the plot. 

Although we have chosen to show only a few repre- 
sentative cases, the bilinearity of all coupling variables is 
well reproduced by our numerical results. It is therefore 
trivial to scale our results to any other desirable values 
of the radial and nonradial average kinetic energies. 



VI. GRAVITATIONAL WAVES 

In order to assess the detectability of gravitational 
waves emitted at combination frequencies, when several 
oscillation modes are present, we focus on CASE IV of 
Table |TI1 where the fundamental radial and nonradial 
modes are present together with their first five overtones. 
The average kinetic energy in the fundamental modes is 
given as m Sec. IVB] with the energy in the overtones 
given again by Eq. and Eq. ([TT]). We also vary the 
average kinetic energy in the radial modes, between a 
minimum value {eI^'°^) = 6.25 x lO"'^ Mq, that resuhs 
in a central energy density variation of 1%, and a max- 
imum value {E^^'^'') = 9.96 x 10^^ Mq, that results in 
a central energy density variation of 5%. We compute 
the characteristic strain of gravitational waves, defined 
as in [H 



he (i^) 



(18) 



where dE / dv is the energy spectrum of the gravitational 
signal and d is the distance of the source. The energy 
spectrum can be written as [47| 



dE 



2)! 



di^ 8 ^ (£- 

e>2 ^ 



2)! 



(19) 



where Zim denotes the Fourier transform of the (total) 
Zerilli-Moncrief function given by Eq. 0. For £ — 2, we 
have 



he (v) 



(20) 



Figure [6] shows hc{i') for a source at the galactic distance 
d — 10 kpc with an initial oscillating configuration given 
by CASE IV and for an emission time of 34 ms. On 
top of this. Fig. [5] also displays the optimal rms noise, 
^ims(^^) for several planned detectors: the interfcromet- 
ric advanced LIGO jd,], the advanced VIRGO .12, 1^, 
the GEO-HF [ij] designs and the proposed wide-band 
dual resonant detector Dual SIC [l^l . We recall that the 
optimal rms noise is defined as 



I (i^) = VShj^, 



(21) 
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FIG. 6: For a galactic source (d = 10 kpc), the figure displays 
the characteristic strain hc{iy) of two initial pulsating config- 
urations for CASE IV of Table HIl In the first {dashed line) 
and second (solid line) configuration, the radial pulsations re- 
spectively lead to 1% and 5% variation of the central energy 
density. The figure also shows the /irms(i^) strain sensitivity, 
Eq. dnj, of AdvLIGO, Adv VIRGO, Dual SIC and GEO-HF 
detectors (see legend). Some of the linear and bilinear har- 
monics are within the sensitivity window of these detectors, 
especially in the Sec ~ 5% case. 



where Sh{v) is the spectral density of the detector strain 
noise. For signals with random direction and random 
orientation with respect to the detector, the rms noise 
should be increased by a factor of \/h [54j . 

Of the several linear nonradial modes excited in the 
simulation for CASE IV, only the two lowest-order ones, 
^f and ^pi have a characteristic amplitude above the rms 
noise of either Advanced LIGO (AdvLIGO) or Advanced 
VIRGO (Adv VIRGO). When a 1% central density vari- 
ation for the radial oscillations is assumed, then all bi- 
linear combination frequencies are also below the strain 
sensitivity, for the above two detectors. However, when 
a 5% central density variation for the radial oscillations 
is assumed, then some bilinear combination frequencies 
are above the rms noise of either AdvVIRGO or both 
AdvLIGO and AdvVIRGO. 

For the GEO-HF detector, a large number of bilin- 
ear harmonics is above the strain sensitivity curve at a 
5% central density variation of radial oscillations. With 
such a detector, even harmonics of up to 10 kHz could 
become detectable. At the lower limit of a 1% central 
density variation, some bilinear harmonics have a char- 
acteristic amplitude which is marginally above the strain 
sensitivity curve. 

Figure [7] focuses on the 1 to 2.5 kHz region of Fig. [H] 
The detectable bilinear combination frequencies are close 
to the frequency of the fundamental ^f mode. Specif- 
ically, the ^P3 — H2 and the ^p4 — H3 combination fre- 




V [kHz] 

FIG. 7: Blow up of Fig.|6]in the region where the chance of a 
detection if higher. Besides the linear ^f and ^pi modes, there 
are several combination tones with a characteristic strain 
larger than the detector noise, in particular for the Sec ~ 5% 
configuration. Note that the combination tones pa — H2 and 
P4 — H3 fall into the same frequency bin and then are indis- 
tinguishable. 



quencies fall into the same frequency bin and their com- 
bined amplitude exceeds the noise curve for both the 
AdvLIGO and AdvVIRGO detectors, for a 5% central 
density variation. It is thus apparent, that some bilin- 
ear combination frequencies may be interesting even for 
the AdvLIGO or AdvVIRGO designs. For the Dual SIC 
detector design, several bilinear harmonics fall within 
its sensitivity window in the range of up to 4.5 kHz. 
Figure [7] shows that the ^P4 — H3 harmonic could be- 
come detectable even at a 1% central density variation. 
Higher signal-to-noise ratios and better event rates can 
be achieved with more advanced designs, such as the de- 
sign for the EURO detector in xylophone mode [s^ . 



VII. DISCUSSION 

Using a gauge-invariant perturbative formalism, we 
studied the bilinear coupling between the radial and non- 
radial polar pulsations in a relativistic star. For typi- 
cal values of mode-energies expected in the post-bounce 
phase of core-collapse supernovae, we find that some bi- 
linear combination frequencies may become detectable 
with planned interferometric detectors, such as Advanced 
VIRGO and GEO-HF, or with proposed wide-band res- 
onant detectors, such as the Dual SIC. 

The possible detection of bilinear combination frequen- 
cies, in addition to the expected linear mode frequen- 
cies, would yield significant new constraints for the high- 
density equation of state of compact stars, as they also 
contain information on the radial modes of the star. Us- 
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ing empirical relations constructed for gravitational-wave 
asteroseismology [sgI , [stI [58| . in combination with nu- 
merical data for radial modes in realistic compact-star 
models [SO^ could allow for the determination of both the 
mass and radius of the star with good accuracy. 

Bilinear combination frequencies of rotating proto- 
neutron stars were recently studied in Q where it was 
found that rotational changes result in the difference in 
frequency between two overtones coming close to the fre- 
quency of a fundamental mode. For such cases it was sug- 
gested that these crossings could result in an enhanced 
gravitational-wave emission at the combination frequen- 
cies, due to 3-mode resonances. The results presented 
here possibly allows for such an interpretation, as the 
amplitude of the ^p3 — H2 and ^p4 — H3 frequencies in 
Fig. [3 is larger than the amplitude of the ^p2 — Hi fre- 
quency. Since the latter involves modes of lower order 
(and of higher average kinetic energy in our examples), 
one would expect the opposite result. It is thus possible 
that the amplitude at the ^pa — H2 and ^p4 — H3 fre- 
quencies is larger because of a resonance effect with the 
linear ■^f mode. Although we have not included back- 
reaction effects, the imprint of a possible resonance may 
be already present in the structure of the equations at 
the bilinear order that is treated here. 

In the present work, we have only focused on one par- 
ticular, cold equilibrium model. In a follow-up study, 
we will survey a larger number of equilibrium modes 
and examine the dependence of the detectability of bi- 
linear harmonics in gravitational waves on the equation 
of state and the compactness of the star. In a more realis- 
tic study, the influence of the high entropy immediately 
after core bounce, the subsequent thermal evolution of 
the equilibrium model and rotational and magnetic-field 
effects should also be taken into account. In addition, nu- 
merical simulations have shown that significant mass ac- 
cretion onto the proto-neutron star (PNS), the settling of 
the PNS to smaller radii (factor of 2) and a non-neglible 
density in the PNS hot envelope, all contributes to the 
gravitational wave signal between the end of the initial 
phase of bounce ring down (~ 10 ms post bounce) and 
the start of convection in the hot bubble (~ 70 ms post 
bounce) [5l| . When additional factors are taken into ac- 
count, our perturbative study could help in resolving spe- 
cific features in the complex gravitational wave signal in 
realistic core-collapse. 

In another study, we are planning on using amplitude 
equations in perturbation theory and a fully nonlinear 
numerical code, as in to also explore possible satura- 
tions of the nonlinear oscillations or the existence of pos- 
sible resonances or parametric instabilities, as suggested 
in 0,&li3. 
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APPENDIX A: NONRADIAL PERTURBATION 
EQUATIONS 

In this section, we write the perturbation equations for 
some of the linear nonradial perturbations of a nonrotat- 
ing spherical star, i.e. the metric variable and the 

velocity 7(0'^) and q;(°'^^. In fact, the closed system of 
perturbation equations for the enthalpy iJ^^'^) and the 
two metric quantities S^O'i) and fc(°'i) is already given 
in go, [HI. 



/(0,1) 



(0,1) 

7,t 



-2$ .e*-^ 



where from one of the TOV equation we have: 



„2A 



(Al) 

(A2) 
(A3) 

(A4) 



APPENDIX B: PULSATIONAL KINETIC 
ENERGY 

The kinetic pulsation energy can be determined as in 
[60| . that in terms of the GSGM variables is given by 



= - / dr {e + p) SuM'e'^ y/^d^x , (Bl) 
2 Jo 

where d'^x — dr dO dcf) and ^ g is the determinant of the 
3-metric. The covariant 3-velocity of the nonradial per- 
turbations is instead given by the following expression: 

Su,, = i i 7(°'i) - ^) eV'", fiC-i) y'™, 0) , (B2) 



where j^^'^') and a^"'^) are two functions that describes 
the radial and latitudinal velocity components. We have 
denoted them with a tilde because they are not gauge 
invariant. However, in the Regge- Wheeler gauge we have 
that 7('''i) = 7(0'i) and = a^°-^\ where and 

are the two polar gauge invariant velocity compo- 



v(0,l) 



nents [32]. When we introduce Eq. (|B2|) into (|BT|) and 
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calculate the integral over the 2-sphere, we obtain 



that reduces to Eq. ([T2)l in the Regge- Wheeler gauge. 



(B3) 
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